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| Abstract 


“ The relationship between booster performance, ion pro- 
 pulsion system characteristics, mission, and payload is dis- 
eussed for the simple case of zero gravitational field. It is 
shown that appreciable increases in payload can be obtained 
_by varying the exhaust velocity in the appropriate way. 


‘Introduction 


_ There has been much discussion in the literature in 
recent years on ion rocket propulsion. Many of the 
problems, such as specific power limitations, thrust 
programming, flight trajectories, and the like have been 
examined quite critically. By now it is well recognized 
that this class of devices possesses inherently low ac- 
celeration but also potentially low mass ratios for 
rather difficult missions. It is not clear, however, how 
such a system might compare with nuclear rockets using 
-hydrogen propellant. To make such a comparison one 
must analyze both concepts in some detail and find the 
areas of superior performance for each. This is much too 
ambitious a task at present since the limiting problem 
areas are not understood in detail for either case. Such 
understanding can arise only after considerable hard- 
ware development has taken place. One can make con- 
siderable headway with respect to the ion rocket, how- 
ever, on the basis of some rather simple ideas. 

We shall consider only the simple case of a field-free 
space. The nature of the specific energy and specific 
power limitations, together with their effect upon the 
relation between payload ‘‘cost”’ and the desired mis- 
sion, will be investigated. The results will illustrate the 
problems involved and serve as a guide for the real 
problems of flight in a gravitational field. 

From the discussion that follows we find that the 
nature of the mission and the characteristics of the pro- 
pulsion plant determine the values for mass ratio and 
exhaust velocity which lead to maximum payload. [x- 
eept for very easy missions (short distances, long mis- 
sion time, or low final velocity) we see that the optimum 
mass ratio is not close to unity. This result is rather 
interesting since the ion rocket is frequently thought of 
as a low mass-ratio rocket because of the potentially 
very high exhaust velocity. 


The Energy Limitation 


Since a nuclear reactor does not at present represent 
an infinite energy source one must consider the specific 


* Work done under the auspices of the U. 8. Atomic 
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energy (total energy available for the exhaust/total 
propulsion machinery mass) limitation and its effect 
upon the payload for a given mission. For very ambi- 
tious missions one may find that performance cannot be 
improved simply by adding more propellant mass, since 
this may imply more nuclear reactor mass to provide 
the extra energy required. Examples will be given later 
which demonstrate that this is a real limit. The energy- 
limited cases correspond to missions for which the 
coasting time is large compared to the powered-flight 
time. In this case, the velocity at the time of cessation 
of thrust is the prime index of performance. 

The problem may be stated briefly: Given a rocket of 
initial mass M, , fuel mass My < M,, and energy E to 
be distributed over M; , which is the maximum value 
of the final velocity V? Let us assume that the exhaust 
velocity may be programmed according to the value of 
m, the total amount of propellant consumed prior to a 
given moment. Then the rocket equation is: 


_ v(m)dm 
oars (1) 
Integrating, 
fy(m) dm 
loan ie oat (2) 


u(m) is an arbitrary function. We want to maximize 
V. — V, subject to the auxiliary condition 


ie MEA aha FE (3) 


The solution of the variational problem posed by (2), 
(3) has been obtained previously.’ 


o(m) = Vago (4) 
where 
M; 
M, = M, — M; 
uw = Mi/M2. 
Substituting (4) into (2) the result is: 
Satine a) 


Compare this with the expression obtained for constant 
exhaust velocity, 


Ve—-Vi=idlnyp (5b) 
Plots of (5a), (5b) are shown in Figure 1. 
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Fig. 1. Dependence of V upon the average velocity 0 and 
the mass ratio u for constant exhaust velocity and for op- 
timum exhaust velocity programming. 


We can now proceed to an analysis of the ‘‘cost”’ per 
unit payload mass and its relation to booster per- 
formance, difficulty of the mission, specific energy of the 
ion rocket power plant, and choice of exhaust velocity 
program. For this purpose, let us assume that the cost 
is proportional to the total system structural mass 
(booster structure plus ion rocket propulsion system). 
This implies that no components are used for more than 
one trip. The major fraction of this cost is then asso- 
ciated with the booster structure (except perhaps for a 
good nuclear heat-exchanger rocket with hydrogen 
propellant) so that the cheapest system results when 
the burnout velocity for the booster is just sufficient to 
place the ion rocket into a satellite orbit. 

The important figure of merit for the booster is then 


B= a My, = initial ion rocket mass 
b (6) 
My = booster dry mass 
For the ion rocket second stage, 


M,=M;+M,+ M,. M; = propellant 


M, 


I 


propulsion CG) 
machinery 


M, = payload 


Then 


ae) ae Y Be GY (8) 
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Define 


My, 
= a ( 


Then, the cost parameter R is 


pa Mot _B8+H1t+3) _ gl 


M L Bd | 
Note that BR is just the ratio of the structural wie 
required for a given mission and payload to the boost 
structural weight required to put the same payload =] 
a low-altitude earth-satellite orbit. 
For an arbitrary exhaust velocity program v(m) we 
may write 


Vem Ve = RG (1) 


h is determined by the choice for v(m). The treatment 
of this section is then more general and may include a 
gravitational field. In the latter case the best choice fon 
the function v(m) would not be Eq. (4). It is con- 
venient to define ‘‘propulsion efficiency” by 


_ (§)M2(V2 — Vi)? _ Ma(V2 — Vi)? 


2 Rae ei 
EB M yw y | bee 1 ( ) 


The justification for this definition is the topic of a 
separate report.’ The specific energy k of the propulsion 
system is . 


k= —, (13) 


where FH is the maximum energy that can be delivered 
to the exhaust gases. The maximum value of k is deter- 
mined by the state of the technology. The parameter 


2k 


sat eee same —1 14 
(V2 = Vi)?’ ? ( ) 


? 
then defines the difficulty of a mission V2 — V,-in 
terms of the energy capacity of the power plant. We 
see that 


M,o e My; 
? = MAVemavon ete. oe 
Using (9), (10), (13), we then have 
1 6 
d= a+8) (16) 


Substituting into (11), we have 


Ke(Bon — 1) = B + yon, (17) 
where Ry is the energy-limited cost index. If we now 


differentiate with respect to » with ¢, 8 = constant and 
set 


we find the minimum value the cost parameter 2 may 


ave for a given ¢, 8. The result is 


(BRe) min. = (18) 


From (17), (18), the corresponding expression for ¢ is 


_at(etn)% 
ou areewe ene ron e LCC 


n* 


d (19) 


Thus, the minimum value R may have for a given ¢, 8 
is given by the parameteric equations (18), (19). 

We can now show that the exhaust velocity program 
(4) gives the absolute minimum for Rz when 8, ¢ are 
given and there is no gravitational field. Suppose we 
have a rocket of mass ratio yw, specific energy k, and 
desire a final velocity V. The energy F is given by 


E=kM,. (20) 


So 
E 
Mi Ml Mis Mai. (21) 


Therefore, if we choose the velocity program which 
minimizes # we have maximized M, and, hence, mini- 
mized M, . This makes R, as small as possible. How- 
ever, the problem of minimizing EH under these condi- 
tions is a variational problem described by Eq. (3) with 
(2) as the auxiliary condition and leads to the same 
solution (4) as before. For the two cases 


v(m) = = V ud (22) 


100 


(B Re) min. 


Fig. 2. Variation of minimum payload cost index with 
nission difficulty for chemical booster (8 « 1). Constant and 
ptimized velocity programs are shown. Energy-limited mis- 
ions. 
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Fig. 3. Variation of minimum payload cost index with 
mission difficulty for optimistic nuclear hydrogen rocket 
boosters (6 = 2). Energy-limited missions. 


and 
v(m) = 5, (23) 


the efficiency functions 7 may be computed from (5a), 
(5b), (12). We have, 
w—l 


7 = (24) 
m 


and 


In*u (25) 


= Lb 1? 
respectively. The parametric equations (18), (19) may 
now be solved. For a given 6, @ (booster performance, 
mission difficulty), (19) gives the value for u which will 
minimize R, . Substitution of this value into (18) then 
gives the corresponding value for R,. Figures 2, 3, 
show the results for 8 = 0, 2. It can be shown that for 
the case v = constant, » < 4.9. This results from the 
fact that » has a maximum at » = 4.9 for this case. 
Higher values of uw are of course possible but would 
lead to increased payload cost. From these figures we 
see that the optimized velocity program gives appre- 
clable benefits in terms of payload cost only for mission 
requirements which approach the energy limits of the 
rocket. 

Now we can put in some numbers to see what con- 
stitutes such an ambitious mission. The most advanced 
current reactor designs achieve about 4000 to 6000 
thermal megawatt days/ton energy release in the core 
before refabrication or other reprocessing of the fuel 
elements becomes necessary. Allowing a factor of 2 for 
the remainder of the structure weight and an overall 
efficiency of 5% for the propulsion plant (this includes 
thermodynamic efficiency of the power plant and ion 
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source efficiency), we have 
k ~ 10" ergs/gm. (26) 


Hence, 
VV <4. 10 em sec: (27) 


This corresponds to missions to the outer planets with 
a requirement that the time be reasonable, say 3 years 
ach way. Technological advances could conceivably 
increase k by a factor of 10 or more. In this case, this 
limitation would be of importance primarily for an 
interstellar mission. 


Specific Power Limitations 


For the 6-year round trip mission to the outer planets, 
specific power limitations may be more important than 
the energy limitations, so we now turn our attention to 
this aspect of the problem. We define the specific power 
K by 

Kas, (28) 

Pp 

where P = exhaust power. The first problem is the inte- 
gration of the equation of motion for arbitrary v(m) 
subject to the condition P = KM,.Weset P = KM, 
since, clearly, maximum acceleration will result when 
the power is maximized. We assume that K is inde- 
pendent of v. If it should turn out that the ion source 
and accelerator structural weights are small compared 
to M/, , then this assumption is reasonable. If the rocket 
thrust is 7’ (a function of the time) then the equation 
of motion is 


dV ii 
= : 2 
dt M,—m ee 
Now 
dm 2p  2KM, 
dt v(m) v(m) 30) 
So 
| v(m) dm = 2KM,t (31) 
0 
gives m(t) for a prescribed function v(m). Also, 
2p 2KM 
4b a ‘p ax Pp bY 
v(m) v(m) ’ iB2) 


which, together with (31), gives LEC s anes (29) 
can now be integrated twice to yield s(r) 
travelled during the burning time 7 hen 


, the distance 


i , 
v(m) dm t 
ae) - i (33) 
2KM, KM, 


Here # is the total energy carried away by the exhaust 
gases and not the total energy capability as in the 
previous section. Therefore / is not fixed. The result of 
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the integration has the form 


PAINE 5 | 
a(r) = eyo oa. ( 
Similarly, 


E _(1+0)2h (+8) 3 
> Rit oR ee Cu Le. 


2 ~ 
s ep 
ie ( ‘ egy) F oe ’ ( 
where F(u) depends only upon the choice of the fune 
tional form of v(m). Clearly, for missions that are no 
energy-limited it is desirable to choose » so that s/ 
is minimized, since this will give the lowest transit tim 
for a given distance s. If the final velocity is to be smal 
(or zero) rather than arbitrary, then the thrust mus 
be reversed at some time 7; , such that V(r) = 0. Th 
form of the result, (36), remains the same, however 
For the choice (4) for v(m) it is interesting to note tha 
dV /dt = constant so that 7, = 7/2. In the followin; 
however, we shall not restrict V(r). 
For a given v(m), the maximization of (36) dete 
mines optimum yu and hence the optimum average ex 
haust velocity 3. Thus, we have a procedure for dete 
mining the best value for 6 when 6 is fixed. Howeve 
the principal problem is the minimization of R, th 
payload cost index, for a specified mission requiremen 
s'/r and specific power capability K. This leads to 
variational problem for the determination of v(m) whie 
has not as yet been solved. The results of the previou 
section do suggest, however, that an interesting clasi 
of functions for the exhaust velocity program is 


100 


(BRp) min. 
S 


| 0 1008 
: 
Fig. 4. Variation of minimum payload cost index with 


mission difficulty and choice of exhaust velocity program for 
chemical boosters. Power-limited missions. 
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8 
Fig. 5. Variation of minimum payload cost index with 
mission difficulty and choice of exhaust velocity program for 
reasonable nuclear hydrogen rocket boosters. Power-limited 
missions. 


v 


(M, — m)"’ Rea) (37) 


v(m) = 


v can be expressed in terms of 7, uw, n, M2 by substi- 
tuting (37) into 


_ KMor 


Pe fers = 
=i IS aie pare (38) 


The resulting form of (36) is 


s eek 2(2n — 1) 
e (1 + 8) ner? — 1) 


ie UG a. | 
E Sher aqy~ i |, 2>0 (39) 


Bie tS 
~ (1 + 8) 


We can now define a parameter 6 which is analogous to 
@ and which determines the difficulty of a mission in 
terms of the specific power capability of the power 
plant. We have 


F(u, 7). 


Gk nie (hts 8) 


40 
oe PF (40) 


6 

We shall see that all missions with the same value for 

6 are equivalent in that they lead to the same optimum 

mass ratio and minimum cost index. The expression for 
the power-limited cost index Rf, is now 


R,(G0F — 1) = 8 + nF. (41) 


Carrying out the minimization process for & as before, 
we obtain a similar pair of parametric equations 


dl 
PF Ue 
(skjia ar (42) 
p/min. dF 
dp 
and 
P+@+u 
a = (43) 


fe? 


For the simple case of constant exhaust velocity 
(m= 0), 


[Gee 2[(u au 5 u] (44) 


Note that the F-functions given by (39), (44) all 
approach the limit (4 — 1)/2 as » approaches 1, as they 
should. The solutions of Eqs. (42), (43), are shown in 
Figs. 4, 5, 6 for n = 0, 1, 2 and 6 = 0, 1, 2. Note from 
Fig. 6 that for v = constant, w < e. This results from 
the fact that (44) has a maximum for » ~ e. No results 
are given for n > 2 since no further performance bene- 
fits result. 

In order to see what these curves mean, let us con- 
sider a specific mission, 


s = 40 million miles 


7 = 65 days 


a 
| 


0.1 kw/kg. 


Then 


| 10 100 


Fia. 6. Variation of minimum payload cost index with 
mission difficulty and choice of exhaust velocity program for 
optimistic nuclear hydrogen rocket boosters. Power-limited 
missions, 
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100 
Ben 
n=l 


LOCUS OF POINTS GIVING 


BRp 
s(km) 
4 Fic. 8. Payload cost index for energy-limited ion rocket 
missions in field-free space. 8 = 2,n = 1. 


In order to obtain these numbers, we find » optimum 
from (43), F from (44) or (39), 6 from (40), M, andl 
M, from (6), (7), (8), (9)- Soptimum is determined by: 
Moptimum, T, K by Eq. (35). The results are shown im 
Table I. The low values for 3 result from ignoring the 
eravitational field. They do indicate, however, that 
for an actual trip from Earth to Mars, the optimum 0 
would be close to AV er for such a trip. Note that mini- 
mizing the payload cost index led to values for 3 much 
less than are usually associated with ion rockets. 


Fig. 7. Nature of the cost index function for a fixed mis- TABLE-I 


sion 0. 
Effect of Exhaust Velocity Program Upon Ion 


Rocket Weight Breakdown 
From Fig. 4 we see that there is nearly a factor of 2 


in payload cost index between the cases n = 0, 2 when Beat x oes ial Seas re 
6 = 0. For 6 = 2 (optimistic nuclear hydrogen rocket se) ae 
booster) the factor is 3. To put the situation differently, 
consider the effect that the choice of n has upon the 3 ; oy Ae a 
mission time for a given cost R, , distance s, and booster yea 0.646 oe 0.287 
capability 6. Take 6 = 2 and assume that 6R, = 10 is M 1.20 1.28 1.30 
reasonable. Then, from Fig. 6, going from n = 0 to M, 0.154 0.355 0.413 
n = 2 reduces 7 by the factor 1.86. Even for BR, = 4 Poptimum » (cm/sec) 8.64 X 10° | 7.96 X 10° | 7.76 X 10° 
the time reduction is 30%. (BRy) min. aa i sa 
In order to illustrate the nature of the minimum in co 7 LAS 
R, for a given 6, BR(w) has been plotted vs p for various Concluding Remarks 


values of 6 in Figure 7 with B = 1, n = 1. Note that 


; ne ; The results of the preceding sections ¢: i . 
there exists a lower limit on yw for a given mission 8, sults of the preceding sestions Gap Desi =a 


even though we have not limited the energy graphically by plotting curves of constant (8R)min, Te- 
: Gh nto UE aged sulting from both the specifi speci 
The 6 Nat (Gor hie Sa a i sulting from both the specific power and specific ener 
We n i‘ a peeehed ings ae nes of conetany ae considerations in the s, 7 plane. This is mie for the 
» note a rather surprising fact: For ‘easy’? missions ral GIL OT OO “Bi : 14 
(0 > 1), the n = 2 case gives a lower cost index than fd corty program. 2 = 04 > onde 
. ‘ 6 * re = ré 2) 1 ‘os res | 
does n = 0 but requires a higher mass ratio p. This ar tae watt/gran in) Wigs (6 and Votes as 
ae ace? i Sof tov ; D Se pe ean distances of the various pl: iN 
results from the fact that for n = 2 the reduction in inciontee oe % We biaing mada a 
fons Sioa eee cm ; sated on the axis of absciss: 
required power plant weight more than offsets the in- Nowe Stee s baci i er 
creased payload and propellant weights. In order to ti Ss a enersy HOE nen only 10h es 
illustrate this point we take 6 = 6.2 and compute pz ai 7 >t > RUS, where ris he ge 
PAN Bono om Ent Cen nate ne the propulsion plant. For ¢ > 7, the rocket would coast. 
y (n and fo1 This situate 
— « ee ieee as bl ‘ S sltue S85 y 1 ‘ie 1 if 
n = 1, 2 from Fig. 6. The booster structural weight is =] ae shown in Fig. 10. Since for m = 1 the 
assumed to be unity. If we further assume + = 65 days aera Se 
ce ee) 
NK 0.1 kw/kg, the corresponding value of @ is also V 
determined. Ve5=2 
5) Qty 
a =Tb 
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s (km) 


- Fie. 9. Payload cost index for power-limited ion rocket 
issions in field-free space. 8 = 2,n = 1. 


he curved part of the (8) min. lines show the coasting 
art of the flight. It should be pointed out that the 
traight sections of these lines are not trajectories, but 
ather the loci of the points s, 7 for which (8R)min. is 
onstant. Since n = | optimizes R for + > 7, while 
= 2 is nearly optimum for 7 < 7, the curves shown 
ere do not give optimum results. However, they do 
show clearly the general nature of the results to be ex- 
pected when an optimization is carried out. 
The picture presented here undergoes major changes 
if real gravitational fields are considered. For example, 
V. — V, is replaced by the characteristic velocity of the 
mission V,. V, is defined by 


Ve [ aw dt 


a(t) = acceleration due to thrust forces alone. 


For values of the acceleration which appear reasonable 
or ion propulsion systems (1 milli-g > a > 0.05 


POWER AND ENERGY LIMITED FLIGHT 
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Frag. 10. Payload cost index for ion rockets in field-free space 


milli-g) V. is independent of a if minimum energy 
transfer orbits are used. In this case, it turns out that 
the equations for energy-limited flight are applicable. 
If, instead, we disregard energy consumption but try 
to minimize transit time, an entirely different sort of 
exhaust velocity program is used. Also, minimizing one- 
way transit time is a different problem than minimizing 
round trip time because of the influence of the motions 
of the planets. In general, then, the situation is con- 
siderably more complicated when actual gravitational 
fields are included. It remains true, however, that ap- 
preciable gains in performance can be achieved by 
proper programming of the exhaust velocity. 
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The Effect of Aerodynamic Forces 
on Satellite Attitude 


D. B. DeBra* 


Abstract 


The effects of torques due to aerodynamic drag and the 
eravity gradient are computed for satellites orbiting at alti- 
tudes of between 80 and 375 miles. Motion about the pitch 
axis is discussed and the equilibrium position is determined as 
a function of altitude for both a dumbbell and cylindrical- 
shaped object. The equilibrium position of the cylindrical 
object in three dimensions is considered as a function of alti- 
tude. Equations are presented in the appendix. 


Introduction 


It is often desirable to place in orbit a satellite which 
is oriented in a unique attitude with respect to the 
earth. If a stable attitude can be found, attitude con- 
trol is reduced to a damping problem. Two methods 
have been proposed to obtain such a stable attitude. 
These methods use, (1) the gradient in the earth’s 
eravitational field, and (2) the aerodynamic torques. 
Kach method requires that the vehicle possess certain 
properties. In gravity gradient stabilization the vehicle 
must have a single axis about which the moment of 
inertia is at a minimum. In the stabilized position this 
axis thus becomes aligned to the vertical [2, 6, 7]. In 
aerodynamic stabilization it is necessary for the center 
of pressure to lie behind the center of mass in the direc- 
tion of vehicle velocity. 


Discussion 


Gravity gradient stabilization can be explained best 
if we visualize a dumbbell shaped object in the earth’s 
field. In a horizontal position the forces on each mass 
(ball) are the same, but in other than this exact position 
there is more attraction on the mass closer to the earth. 
Therefore, in addition to the net force there is a moment 
acting unless the dumbbell is aligned to a vertical 
position. 

It is difficult to design a vehicle without consideration 
of both gravity gradient and aerodynamic torques. 
When attitude stabilization is attempted utilizing 
either effect, disturbances may be introduced from the 
other. The problem considered here is a study of the 
satellite stable equilibrium position in a circular orbit, 
under the influence of a gravitational field, moving in a 
model atmosphere. 

Two nominal vehicles have been chosen to illustrate 
this problem and are shown in Fig. 1 and 7. The ge- 
ometries are distinct but for simplification the physical 
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properties chosen are the same. The analysis procee 
from a consideration of behavior about a single axis, thi 
pitch axis (discussed in ref. 1), to a three dimensional 
study of the cylinder. 

Aerodynamic forces at altitudes above 80 statute 
miles make an interesting study in themselves. Th 
computation of drag based on free molecular flow i 
complex [3]. The most exact theory takes into account 
not only the random velocity of the particles as they 
approach the vehicle but also as they leave the surface 
after impact. The molecular speed ratio is the ratio of 
vehicle velocity to the mean particle velocity. When 
this ratio is high, the effect of the incoming particles 
can be computed quite accurately by ignoring the 
random components. When random motion is only con- 
sidered as the particles leave the vehicle surface, the 
theory is described as Newtonian-diffuse. 

Particle random velocity defines the ambient tem- 
perature. At these altitudes the temperature may be 
around 1500°R. After a particle strikes a surface it 
loses some of its energy and leaves with a velocity as- 
sociated with a temperature near the surface tempera- 
ture of the vehicle. This may be around 550°R. It 
would seem natural, then, to neglect the momentum 
change of particles leaving the surface since their mean 
velocity is even smaller than the random velocity 
ignored in the Newtonian-diffuse analysis. But, in some 
cases, this is not true since the particles’ effect is as- 
sumed normal to the vehicle surface. The relative force 
may be small but the resultant moment can be signifi- 
cant. When the effect of particles leaving the surface is 
also ignored the flow is Newtonian and the computation 
is considerably simplified. 
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FIGURE 2 


Newtonian flow possesses a drag coefficient of two. 
hus, when the drag coefficient based on free molecular 
ow approximates two, the assumption of Newtonian 
ow yields reasonable results. Figure 2 presents a 
nethod of determining in what range of altitude this 
pproximation can be made. The drag characteristic 
slotted is for a sphere but the shape of the curve is 
haracteristic. The assumption of Newtonian flow is 
herefore justifiable for altitudes less than 375 statute 
iles for accuracies near 10 percent. 

The Knudsen number is the ratio of particle mean- 

ree-path to a significant vehicle dimension. Irom ex- 
berimental results flow is essentially free-molecular for 
udsen numbers greater than two. For the vehicles 
nder consideration here, the Knudsen number is 
pproximately five at an altitude of 80 statute miles 
nd increases rapidly above that. 

_ All computations in this paper are based on New- 
onian flow and are a good approximation for altitudes 
etween 80 to 375 statute miles. 

First let us consider the vehicles of Tig. 1. With 
reedom of motion only about an axis normal to the 
lane of the orbit, the variable describing the attitude 
is merely the angle a defined in Fig. 1. 

Consider a vehicle stabilized by the gravity gradient 
so that its position corresponds to a = 90 degrees. lor 
the assumed separation of the center of mass and the 
center of pressure of 0.1 ft, the torque required of a con- 
trol system is shown in Fig. 3. This assumes that a 
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constant torque must be produced by the control 
system in order to maintain position. If not, the vehicle 
will rotate toward the position in which the gravity 
gradient torque balances the aerodynamic torque. 
Since a control system that produces a constant torque 
is expensive in weight and power, the attitude error 
that can be tolerated determines the minimum altitude 
of a circular orbit when constant control is not used. 

For the model considered, it may be correctly argued 
that the deviation from the vertical can be planned on 
and should not be considered an error. For a satellite, 
however, the aerodynamic forces would not be as con- 
stant as they are in the model. Normal fluctuations in 
the atmosphere and variations in the vehicle altitude 
cause a variation in the equilibrium position. Since the 
variations may be close to the deviation from the verti- 
cal, where a nears 90° the deviation is treated as an 
error. 

The aerodynamic pressure decreases with altitude 
but theoretically never vanishes. Therefore, if there is 
an area moment about the center of mass that can pro- 
duce an aerodynamic torque, there is no altitude at 
which the equilibrium attitude of the vehicle exactly 
coincides with the vertical. However, the error becomes 
insignificant for higher altitudes. 

When the vehicle is stabilized by aerodynamic effects 
(a = 0) at a given altitude, the restraining force is a 
function of the deviation. The frequency of oscillation 
about the equilibrium position is a convenient method 
of comparing the stiffness of constraint for small devi- 
ations. 
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Vigure 4 presents the natural frequencies of these 
vehicles for small oscillation about equilibrium posi- 
tions determined in the absence of all torques except 
the constraining torque indicated. 

It is characteristic of the gravity gradient torque to 
vary approximately as the sine (2a). In other words, 
there is no constant disturbing torque due to the gravity 
gradient when one of the vehicles is in a horizontal po- 
sition (a2 = 0 degree or 180 degrees). This does not 
mean that the position is stable. When a is near 180 
degrees, the slightest disturbance allows both the 
aerodynamic and gravity torques to increase the devi- 
ation. This condition is clearly unstable but the position 
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near a = O degree is not. At a near 0 degree, the re- 
straining torque due to aerodynamic forces may be less 
for small deviations than the gravity gradient divergent 
torques. 

The crossover point (where the aerodynamic restrain- 
ing torques equal the gravity gradient disturbing 
torques) is an altitude below which the position, a = 0 
degree, is a stable equilibrium position for aerodynamic 
stabilization. By the nature of the gravity gradient 
torques the effect near a = 0 degree are equal in magni- 
tude to those near a = 90 degrees, but reversed in sign. 
Therefore, Fig. 4, which displays the natural fre- 
quencies, also indicates the altitude at which this cross- 
over occurs. 

For the vehicles chosen in this part of the study, 
aerodynamic stabilization can be accomplished without 
theoretical error for orbits below 120 statute miles. 
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The two important considerations for a stabilized 
vehicle are its equilibrium position and how stiffly it is 
held in this position. A method of displaying both is 
shown in Fig. 5, where the net torque acting on the 
vehicle is plotted as a function of vehicle attitude for a 
certain altitude. The stable equilibrium position is the 
point where the net torque passes through zero with 
positive slope, and the slope indicates the stiffness of 
constraint. 

Figure 6 shows the equilibrium position as a function 
of altitude. The rapid change in position with altitudes 
near 150 statute miles should be noted. Variation in the 
atmosphere and in vehicle altitude would have the 
greatest effect on the equilibrium position at this 
altitude. 

An additional source of torque must be considered 
when the three dimensional case is analysed. It is a 
gyroscopic or precessional torque. The reference frame 
used in the previous analysis involves the vertical which 
rotates in space as the satellite circles the earth. A body 
which is fixed with respect to that reference has the 
same angular velocity. Components of the angular 
velocity along the principal axes define angular mo- 
menta, which, in general, are changing direction. For a 
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given position, this results in torques about the body 
axes and a coupling of the motion between the axes. — 

One effect is to make the equilibrium position of th 
satellite unique when the moments of inertia about t 
principal axes are distinct. This can be verified by 
simple experiment. If an object with three distin 
principal moments of inertia, such as a book, is spu | 
about a principal axis as it is thrown into the air, it is 
found that the motion is stable about the axis of 
maximum or minimum moment of inertia, but unstable 
about the axis of the intermediate moment of inertia. 
Apply this to a satellite in the absence of aerodynamic 
torques when the satellite axis of minimum moment of 
inertia is constrained to the vertical by the gravity 
gradient torque. Then, to attain attitude stability, the 
axis of maximum moment of inertia must be aligned 
with the orbital angular velocity vector. 

To study the effect of aerodynamic torques when 
there are three degrees of angular freedom for the 
cylinder, consider the previous analysis to be altered by 
changes in the cylinder, as shown in Fig. 7. The change 
in the moments of inertia reduces the gravity gradient 
torques by a factor of nearly two and the location change 
of the center of pressure produces aero torques about 
the other axes. The resultant equations, which are 
highly non-linear, are developed in the appendix (10). 

There are many ways of defining a vehicle attitude. 
The method used here is common and involves the 
angles as defined in Fig. 7. The effects of the atmospheric 
torques are described by the behavior of these angles as 
a function of altitude, as shown in Fig. 8. 

There are two results of primary interest. One occur; 
when the aerodynamic effects have tipped the vehiel 
over to a horizontal position and the axis of minimun 
moment of inertia is no longer aligned to the vertical 
The satellite does the “next best thing” and aligns the 
axis of intermediate moment of inertia to the vertical 
The second is that for small deviations from the aero 
torque-free position, the equilibrium position ean bi 
computed roughly by using the spring stiffness giv 
for each axis (11), : : ee 
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Conclusions 


It should be emphasized that the results presented 
here apply to certain vehicles and differ greatly as the 
vehicle characteristics change. For instance, the ratio 
of gravity-gradient torques to aerodynamic torques 
varies as the square of the size for similar vehicles, and 
greater effects can be expected if specific properties of a 
vehicle are changed. 

If a vehicle is designed specifically for gravity 
gradient stabilization, it is likely that a center-of-mass— 
center-of-pressure separation could be controlled to less 
than two percent of the vehicle length. Similarly, aero- 
dynamic stabilization could be achieved more easily in 
a vehicle which is designed with nearly constant 
moments of inertia about all axes. 

Although any position of the satellite vehicle may be 
maintained by the use of an attitude control system, the 
power requirements for such a system could be pro- 
hibitively high in any satellite designed for long orbital 
life. For vehicles designed for a stable equilibrium po- 
sition—so that the potentially complex control function 
is reduced to a damping problem—allowances must be 
made for the effects of extraneous torques which could 
change the equilibrium position beyond the design 
tolerances of the vehicle attitude. 
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Appendix 
Nomenclature 

L Angular momentum of satellite 

N Torque acting on satellite 

W Angular velocity of satellite 

h Displacement of center of pressure re. the 
center of mass 

F Force 

a Angle of attack of the body axis of symmetry 

[A] Matrix defining transformation from orbit 
coordinates to body coordinates 

Qj Elements of [A] 

[A ;5] Matrix of transformation of jth triad into 
(7 + 1)st triad through an angle A;; 
about the 7th component of the jth co- 
ordinate triad. 

Ci; Cosine of angle A ,; 

Si; Sine of angle A;; 

L six The 2th elements of the intertia tensor 
expressed in the kth coordinate triad 

I; The moment of inertia about the 7th 


principal axis of the body. 

= Denotes a vector 

ij A unit vector along the 7th component of 
the jth coordinate triad. 


The coordinate triads are defined by the subscripts 0, 
1, 2, b which define coordinate triads successively from 
the orbit triad after transformation [A 101, [An] and 
[A32] as shown in Fig. 7. The final coordinate triad, 6, 
corresponds to the principal body axes in which the 
magnitude of the principal moments of inertia, 1; , Lo, 
andl, (lei, <J,). 

The components of triads are numbered successively 
1, 2, 3 and form a right-handed coordinate system. 

Subscripts a and g refer to aerodynamic and gravita- 
tional effects. 
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Equations of the Applied Torques 


The basic equation for angular motion is 


where 


d d os 
+ Ref.9) (2 
(Sue Cet +W-X | (Ref.:9)aa 


Since the equilibrium positions are sought the trans- 
formation [A] is constant and 


d 4 
lees. = 0 (3) 


in both body and orbit coordinate frames. Computation 
in principal body axis results in considerable simplifica- 
tion because the inertia tensor is diagonalized and 
moments of inertia are constant. The two sources of 
torque considered are (1) Aero dynamic torque and 
(2) gravity gradient torque acting on the cylinder of 
Eigva 

(1) Aerodynamic Torque (For Newtonian diffuse 
method, see [3]) 

Assuming Newtonian flow, a center of pressure and 
therefore an h exists for the vehicle of Fig. 7. Then 


Fy = Fx20 = —pS20 (4) 
and 
0 dye | 
FP, = [A] | Fo} = — pS | dg }. (5) 
0 a | 
Since cosa = ay, and S = reosa + 10sin a, the 


ojected area can be found and p, the aerodynamic 
essure, is determined from Fig. 3. Therefore 


hsd22 
(N53 = (hy Xx faoys = pS| hay — hsay (6) 
—hyar2 


(2) Gravity Gradient Torque 
This torque is computed from the potential energy 
e body possesses as a result of its attitude, i.e., the 
otential can be divided into a potential of the total 
ass if located at the center-of-mass plus a potential of 
e body about its center. A spherical mother body is 
ssumed here but the method is the same when oblate- 
ess is included. This has been done in [7] and is briefly 
utlined below. 

To find the total potential for a rigid body, the unit 
otential can be expanded in a Taylor series and in- 
egrated over the body. The significant term of the 
otential using Wo = go/Ry and coordinates of Fig. 9 is 


— (Ie + Is30 — 2I110)Wo /2 


The introduction of the moments of inertia allows a 
nvenient method of expressing this part of the po- 
ntial as a function of the angles of a coordinate trans- 
rmation. [8 and 9} 


[To = [AULA] 


~The gradient of a potential is a torque when the co- 
rdinates are angles. The component torques are about 
he axis for which the angle rotation is defined. 

It is appropriate here to discuss the choice of co- 
rdinates to define a satellite attitude, that is, to choose 
he form of the transformation [A]. Two types are 
ommon: (1) Classical Kuler angles are defined by three 
otations, the last of which is about the same coordi- 
ate axis as the first [9]. There are many variations in 
he choice of first and second axis, but all have the 
bove property; (2) The only distinct procedure for a 
hree-angle transformation is to have one rotation 
bout each of the coordinate axes [8]. 

For either type of transformation it is most con- 
enient in this application to choose the first rotation 
bout the axis of symmetry of the potential function, 
he vertical. Because of the symmetry, the potential 
imction will not be a function of rotations about this 
xis; therefore the labor in finding the torques is re- 
uced. 

When small motions are instructive in understanding 
ne physics it is easier to use the second transformation. 


Therefore defining 
[A] = [Ax|[Aa][ Aro] 


the torque in terms of the principal moments of inertia 
is 


SoCoSio(L1C32 + Iz S32 — Is) 
(N,), = 3Wo | SuCuCm(iCie + IxS32 — Is) | (7) 
Sx2Cs2C21(L1 — Iz) 
Equations of the Equilibrium Position 
Substituting (2) and (3) into (1) 
WXL=N,+N, (8) 
letting Ly, = Wy ete. and 


Ws = [A|Wo = Wo Ea (9) 


then (8) becomes 
WaWw ls — Iz) + (Na)w + (No)w = 0 
WaWrUs — li) + (Na)w + (No) = 
WwW — Ie) + (Na)e + (No) = 0 


| 
co) 


(10) 


where the components of V,, NV, and W are given in 
(6), (7), and (9) as functions of Ay, Aa, and Az . 


Small Angle Expressions 


The physical cross-coupling due to the angular 
velocities relative to the orbit coordinates vanishes for 
an equilibrium position when the vehicle is at rest rela- 
tive to the orbit coordinates. The coupling with the 
orbit rate (due to the rotation of the orbit coordinates ) 
still exists. This coupling combines with the gravity 
gradient to form an effective spring restraint from 
which small error angles due to external torques can be 
estimated. For the aerodynamic torques the error angles 
are approximately : 
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First Order Error Propagation ina 
Stagewise Smoothing Procedure 
for Satellite Observations | 


Peter Swerling* 


Abstract 


A practical method of smoothing satellite data by evaluat- 
ing a finite number of parameters, or elements, is presented. 


Introduction 


The subject discussed is that of smoothing observa- 
tional data in cases where the observations, in the ab- 
sence of observational error, would all be determined 
by the time of observation plus a finite number of 
parameters, called elements. The objective is to esti- 
mate the elements. 

The immediate motivation for this arises in studies 
of estimation of earth-satellite orbits from observa- 
tional data. In this case, if, for example, the force field 
were known exactly, one could regard the elements as 
the position and velocity components at a particular 
instant ¢ ; if the field were a central inverse square 
field, the elements could alternatively be the six ele- 
ments of a Keplerian elliptic orbit. (It should be men- 
tioned, however, that satellite orbit prediction is only 
one of the possible applications of the results. ) 

In satellite tracking and prediction, it is desired to 
produce ephemerides—i.e., predictions of the future 
position as a function of time—as well as to make 
various other types of decisions and predictions. As 
new observations become available, one can improve 
the accuracy of these predictions. Ideally, one would 
like to store all previous observations of the object, 
and combine these in some optimum way to yield the 
desired predictions or decisions. The optimum method 
for processing the available data would be based on 
analysis of the error statistics for the individual ob- 
servations, and on the functional dependence of future 
predicted quantities on the previous observations. 

In satellite tracking one is dealing with a situation 
in which there may be a large number of observations 
of varying degrees of accuracy, as well as large numbers 
of tracked objects. Also, methods of orbit prediction 
(even in the absence of observation errors) are subject 
to various sources of error, such as 

(a) uncertainty in the forces acting on the body 

(earth’s gravitational and magnetic fields, at- 
mospheric resistance) ; 

(b) cumulative errors in solving the equations of 

motion. 
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Two features would be desirable in a tracking and 
prediction method: - 

(1) The data processing load per tracked object. 
should not exceed a certain maximum, regardless of 
how many observations are available to be processed. | 
On the other hand, the prediction method should not 
throw available observations away. | 

(2) The method should be adaptable to situations 
in which the underlying prediction functions are subject 
to the above-mentioned uncertainties. 

The stagewise procedures described below are moti- 
vated by these considerations. 

The particular methods of data smoothing to be dis- 
cussed are variations of the classical method of minimiz- 
ing a quadratic form in the residuals (in practise this is’ 
usually a weighted sum of squares of residuals). After 
defining this method, the first order dependence of the 
errors in the resulting element estimates on the ob- 
servation errors is established. 

We then go on to describe a stagewise procedure for 
processing the observational data, in which the element 
estimates at each stage are smoothed in a particular 
way with some additional observations. This is, in es- 
sence, a type of differential correction procedure. It is 
shown that the errors in the resulting element estimates, 
after stagewise smoothing of a given set of observations, 
have the same first order dependence on the observa- 
tion errors as would the errors in the estimates obtained 
by simultaneous processing of the same total set of ob- 
servations. 

Some statistical properties of errors in the element 
estimates are derived for the case in which the observa- 
tion errors are regarded as statistical variables and in 
which the matrix of the quadratic form to be minimized 
is the inverse covariance matrix of the observation 
errors. (For Gaussian error statistics, this would result 
in & maximum-likelihood method of estimation. ) 

The elements are at first regarded as constants, and 
later the treatment is extended to the case in which the 
elements are regarded as time-dependent. 

We assume there are N observed quantities FF", 


w= 1,---, N, each F“ being a real scalar. We also as- 
sume that n real constants x;,7 = 1, --: , n exist such 


that in the absence of observation error, all N observed 
quantities would satisfy the relations 


P= Eas > °°" Un, tu) (1) 


The vector x = (2, -++ , a) will be called the “element 
vector,” and its components the “elements”; t, is the 
time at which the ut” observation is taken. 

When observation errors are present, the 7 are 
given by 


eat (ppm arn eiey a ee. (2) 


For purposes of illustration, consider the observation 
of a satellite following a Keplerian orbit. The elements 
might then be taken to be the eccentricity, semi-major 
axis, inclination, etc.; the quantities F“ might be ob- 
servations of such ae as range, azimuth, elevation, 
or range rate from particular observation sites; the f” 
would be the functions describing the dependence of 
these observed quantities on the elements and time; 
and the &" would be the observation errors. 

The times ¢, need not be in any particular temporal 
order, nor are they necessarily all distinct. They are 
assumed to be measured without error in all cases in 
which the functions f“ depend explicitly on ¢,. There 
is no loss of generality in this assumption, since timing 
errors may always be reduced to equivalent observation 
errors; this would of course modify the statistics of the 
resulting equivalent observation errors (for example, 
this might introduce an additional systematic com- 
ponent into the equivalent observation errors). 

_ Many of the formulas below will involve the functions 

f* and their partial derivatives df"/dx,;. Practical ap- 
plication of these formulas would be possible both for 
cases where analytic expressions are known for the f* 
and their partial derivatives and for cases where these 
functions must be evaluated by numerical integration, 
as well as for cases where some of the functions have 
known analytic expressions and others must be de- 
termined by numerical methods. 

Henceforth it will be supposed that N 2 n. The 
problem to be considered is the estimation of the ele- 
ments by means of smoothing, in some sense, of the ob- 
servations. 

A classical smoothing method is as follows: writing 
= (P,,:--, Pp) for an estimate of the element 
vector, and writing f(a, t) for f(a, ---,%,t), f(P, t) 
for f(P:,---, Pn, t), and so forth, the method con- 
sists of minimizing with respect to P:,---, Pn the 
quadratic form 


N 
= De ul EY 
Bjv=1 


where (7) iS a symmetric, positive definite matrix. 
Thus, the method consists in minimizing a positive 
definite quadratic form in the residuals. Differentiating 
Q with respect to P; and setting the results equal to 
zero, we find that the minimizing estimates P; must 
satisfy 


= f'(P, i Ie — ~-, t)] (3) 


— Of” : a 
ares aye CP, t, ) [F" fo (RP; ta) ln = 0 (4) 


(¢ = 1; <=: 7) 


It is clear that if &” = 0, all uw, then P = x is a solu- 


tion of (4). The first question to be investigated is that 
of the first order propagation of errors—i.e., the first 
order dependence of P — x on the errors &". It will be 
assumed that the functions f” are sufficiently well be- 
haved for the following operations to be valid. We may 
write 


PT saa Orne) 


+pe CHMly Ee eee” 
© (P,4) = © (4) 
; (6) 
+ oF waa) + 


Neglecting all terms of higher order in P — x, and 
substituting (5) and (6) into (4), we obtain 


2 cal (a, 6) as af’ ee tes -x)| 


B,y=1 ji Ox; , OV; 
(7) 


x| P - #*t2, 4) - Ee (x, ty) (P53 — ») |=0 


It is also clear that for sufficiently small &" and 
P; — x;, the term involving second derivatives of f” 
may be neglected. The result may conveniently be ex- 
pressed as follows: let 


r) Bu 
ata, 4) = © (x, 4) (8) 
OX, 
N 
pi(x) = ys wi (x, ty) [FY — f(x, ty] (9) 
Lv 
p(x) = {pi(x), ee. pn()} (10) 
N 
By(x) = oe, thins CL, by )Gy (ly dy) (11) 
bv 
Bw) = {B;;(x)} (ale) 
Then, assuming B(x) to be non-singular, 
P — «x = [B(x)J'p(2) (13) 
Eq. (13) expresses the first order dependence of 
P — x on the errors &". This can be expressed equiv- 
lently: 
N 
Pr 2, = Dlr) FS fe (14) 
B= 


Ti(2) = 3) DY (BG) inal, 6) (5) 


A special case of this is as follows: suppose one has 
already an estimated element vector p = (1, °** 5 Dn) 
together with K new observations. One can form a new 
estimate vector P by smoothing the original estimate 
vector p with the new observations in the following 
manner: P is determined by minimizing, with respect 


to P,,-+--, P,, the quadratic form 
= tur Pu — Py) (p» — P,) 
‘ (16) 
n-+K 
+ DU tol — FCP, ta) LP” — FCP, te) 
a) 
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In this case, 


f(z, ti) = tu, w=, ey 
af (a, ty) = Ou le a (17) 
Ms, ils ‘ »n 
ithe = Pu b= 1, i ,n 
Also 
pi(x) = > Pg Oe 29) 
rs 
ae (18) 
ae MS Nut” = cise tule (2, ty) 
Byv=n+1 
(4 = ike ? i) 
and 
n-+K 
Die Te 
Bi(@) = ag + -s uv (2, ty )as (2, ty) (19) 


(a7 = its hg , 1) 
If we also define r"(p, t,) foru=n+1,---,n+K 
by 
(Dp, ty) = BS — fp, te) 


(20) 
pant il,->-,n + Kh 
then p;(2) becomes, to first order, 
p(x) = Dy Bijs(u) (pi — 25) 
j=1 21) 
n-+K (2 
aa De Nw; (2, Ga, hi 
B,y=n+1 
Consequently, (13) reduces to 
P—z=p—2£+ [B(x)] p*(x) (22) 
where 
n+-K 
pi (xv) = Der Nuvhy (0, tr Cp, ty) (23) 
Lyn 


Hq. (22) may be used as the basis for a first order 
differential correction to p, given the additional observa- 


tions F", 4 =n+1,---,n+K. First rewrite (22) as 
P— p= |B@)\"p*(@) (24) 

Then, to first order, one can also write 
Ep (Bip) pp) (25) 


Since all the quantities on the right of (25) are known, 
(25) gives the required first order correction to p. 


A Stagewise Smoothing Procedure 


Suppose the matrix (7,,) can be written as a diagonal 
array of matrices 


(1) 
U} 


| awe (26) 


(s) 
ul 
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where 7, s = 1,--:, S, are Ns X N,; matrices. Fc 
future notational convenience, we ill regard the in- 

& 8 ay 
dices of the components 7, of 7 as running Ove 


p, v= Mis, a Lee , M, where 


DN, 


i r=1 

M, = 0 cys) 
M, ze ete, 

Mg = N 


We define a stagewise smoothing Dean as follows 
Suppose initial element estimates } } are obtamg 
by minimization with respect to P of the quadratie 
form 


; i 
g® = 2 nek =f (Poel (28) 


(Ff CP ae 


Now define sequences of matrices B® and B® (aya 
follows :* 


(s-) 
BY = Bi; 


BY (2) = BS? (x) 


Ms 
4 Pate aia eas 
4.v=M,—1+1 | 
| 
BS =0 (29) 
} 
: 


Ms 1 
ar SG ns a; (2, BYES Ge tu) 
H,v=M s_j+1 


BE (2) = 0 


k : s) a ; 
For s > 1, the s"" element estimates {P\} are ob- 
. Jee a Oars ‘ (s) 
tained by minimizing, with respect to P™’, the quad- 
ratic form 


n 
Om = S. BY PS = Per id P| 


t,j=1 


+ Do gree, a 


B,Y=M 54+ 


eae A Cada” 


Let us also define p$ (x), s = 1, 2, --- , S, by 
pi (nie Do BES GP a ted 
j=1 
ae (31) 
— »s nae a; (x, oe Oli <F Pe tu) | 
B,Y=M ._4+1 


Then, using (13), (18), and (19), it can be shown that 
to first order 


P® i [B® (2) Tp (x) (32) 


Alternative definitions for the matrices B® are possible 
see the remark at the end of this Section. 


It is also not hard to show by an vinductive argument 
that 


PO — x = [B® (x)J "9 (2) = [B(a)Po(z) (33) 


where B(x) and p(x) are defined as in (9), (10), Gul) 
(12). This is proved by showing that B® (x) 
e(x), and p(x) = p(z). 

Thus, the first order dependence of the errors in the 
estimates {P{°} is the same as that for the {P} ob- 
tained by processing all N observations at once by 
minimizing @ as defined by (3). Another way of stating 
this is to say that, to first order, the estimate obtained 
by processing the N observations by the stagewise 
procedure just described is the same as that obtained 
by minimizing Q as defined by (3). (However, the 
range of magnitudes of 8" for which the first order 
expressions give good approximations to the estimation 
errors is not necessarily the same for the stagewise 
method as for the non-stagewise method. ) 
_ A stagewise smoothing procedure may be advan- 
_tageous in certain situations. For example, suppose 
that observations are coming in at some average rate; 
in the stagewise procedure, it is not necessary to store 
all previous observations and N xX N matrices (m,) 
with N increasing. It is at any time necessary to store 
only the ‘current’ element estimates {P{°} and the 


matrix B, that is, 5 (n + 3) quantities (taking into 


account the symmetry of B). 

It can be seen that the matrices B® play the role of 
estimates of the matrices B® (a). Thus, the particular 
method of defining the sequence {B“} above is not 
the only one possible. For example, one could define 


BS = BY, ¢=1 


Ms 
) 
vase pee 


H,v=M 5_1+1 


Gees, ty) 


igor tl id Bie aay (29') 


se dl 


and define the stagewise process using the matrices 
B® instead of B®. Comparing with (29), we see that 
the main difference is that B® can be computed, for 
s > 1, before one computes Pp”, 

This lends itself conveniently to first order determi- 
nation, for s > 1, of P® — P°” by means of Kq’s. 
£20),. (23), and (24). One would write p = P“”; 
(ae P®. and Pp” — per 1) [Be 1 Be ( pe “i 
(The vector p* would be defined by obvious modifi- 
cations of (20) and (23).) 


Application when the Elements are Functions 
of Time 


Suppose the elements x, ---, 2» are functions of 
time: x; = 2;(t). Also suppose that error-free observa- 
tions are given by 


P= g{x(tu), tl (34) 


Suppose the manner in which x(t) depends on its 
values at some t is known: 
w(t) Pa F[x(to), to 5) t (35) 


Then we may obtain the case considered in the pre- 
vious sections by defining the elements x; in the for- 
mulas of those sections to be x; = x;(t)) and by defining 


Tes tu| = g'{Filx, to ) WAL at ase) Slay to , bls tu} (36) 
(where « = x(t)). 
Modification when the Functions f" are [mper- 
fectly Known 


Suppose the observations are given by 
PY = hx, ty] + € (37) 
(where &" are observation errors), 


but that the estimated elements {P,} are obtained by 
minimizing Q as defined by (3), with the functions f” 
(which differ from h*). 

So long as F“ — f* are sufficiently small, equations 
(9)—(18) still describe the dependence of P — x on 
Fr" — f". This dependence was expressed 


| Peon ee p> jae a pinee tu) | (14) 


where: 
N 


Beene, ij Mav j aes ua) (15) 


Per) = > 


The observation errors &” are now given by 


i a (38) 


Vv 


Therefore 
"f(a, ty) = 8 + Wa, te) — fC, ty) (39) 
Thus, (14) becomes 


YT (a)e 
+ pai 


Statistics of Propagated Errors 


(x) [h"(x, ty) — fu(z, tu)). (40) 


Eq. (13) may be used in an obvious manner to de- 


termine the means and covariance matrix of {P; — x} 
as functions of the means and covariance matrix of 
{e}: 


We shall deal here with a special case, namely, one 
in which the ensemble means and covariance matrix 
of {8"} are known and in which the matrix {7} of the 
quadratic form Q has a special relation to the covariance 
matrix of {&"}. 

Since the ensemble means of {8"} are assumed known, 
it is no loss of generality to assume they are zero. In 
this case the covariance matrix is (denoting ensemble 
means by E(_ )) 


du = H(8'8’) (41) 


The Journal of the Astronautical Sciences 49 


It will now be assumed that the matrix (mu) in 
(3) is 
(mu) = (¢y) | (matrix inverse ) (42) 


If {8"} were to have a Gaussian probability distribu- 
tion (and if f“ = h”, all w) then the resulting method of 
obtaining P would be the maximum likelihood method. 

The covariance matrix of {P; — x} assumes a par- 
ticularly simple form in this case. For generality, we 
will deal with the case described in the last section, in 
which f“, the functions used in the quadratic form Q, 
may differ from h". 

The means of {P; — x;} are 


KP; — a) = TFA, ) —P@, 6) 43) 


The covariance matrix of {P; — aj} is readily estab- 
lished to be 


EXP; Cae A ae E(P; rad xi)} 
-{P; — 2; — E(P; — 2;)} = [(B@)li 


This holds whether the observations are processed 
all together or by a stagewise procedure as described. 
In the latter case, of course, it must be assumed that 
the covariance matrix of {&"} can be broken up into a 
diagonal array of covariance matrices corresponding 
to the different stages. It is, in fact, quite easy to 
establish that B® (a) is the inverse covariance matrix 
Of (Pe a}: 


Be te Ce ee) 


45 
-{P;? — 2; — E(P;? — 2;)} = (BY (ala au 


This throws some further light on the stagewise 
method of Section II. The matrices B“” occurring 
in the quadratic forms Q™ are seen to be estimates 
of the inverse covariance matrices B“ (x) of the 
element estimates P“”. Thus, if the error statistics 
were Gaussian, this procedure would consist at each 
stage of a maximum likelihood smoothing of the pre- 
vious element estimates with the new observations. 

The above formulas may be used to determine the 
rate at which the covariance matrix of the errors 
P{? — 2; decreases as additional observations are 
processed. In fact, this information is contained in (44) 
and (45). 

As a special case, consider the case where all ob- 
servation errors are mutually uncorrelated. (If this 
is not true originally, it can be made true by means of 
linear transformations.) In this case, the matrices o 
and 7 are diagonal. 

If we now regard P“’ as the element estimate vector 
resulting from the processing of the first s observations, 
we have the inverse covariance matrix for ee ae x; 
given by 


Bi? (x) = p> Nes (ie, tu On Meet (46) 
= 
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(This holds whether the s observations were processec 


all together or stagewise. ) ; 
It is also quite easy to verify that, for s 2 n, 


[B® (x)|g = (Bo? @) li — d(x) d(z) (47 
where 
ds” (x) ; 
/ qu Dd, Cet IBoe en 
7 4/3 2 Nss Ss [Bo (a) ]7¢ anf (2, sy) ae, ts) 


ik=1 


(48 


We might close this section by briefly discussing th: 
subject of systematic errors, which can be defined as 
errors which are highly correlated over a certain group” 
of observations. There are several possible ways of 
handling systematic errors. They could be regarded 
as introducing off-diagonal elements into the correla- 
tion matrix ¢, and the smoothing matrix 7, , the 
magnitude of the off-diagonal elements being chosen” 
according to the approximate magnitude of the syste- 
matic errors. Or, alternatively, a systematic error 
might in some cases profitably be regarded as an 
additional parameter to be estimated from the data— 
i.e., as an additional element. This is one of the choices 
that would have to be made in actual implementation 
of these results. 


Modified Stagewise Procedure for Time-Varying 
Elements 


In this section the elements will be considered func- 
tions of time, x; = 2x;(t). It will be assumed that the 
elements vary much more slowly with respect to time 
than do the functions describing observations. 

We suppose that the matrix (7,,) can be written (as 
in Eq. (26)) as a diagonal array of matrices 
n,s = 1,2,--- ; the further assumption is made that 
the times ¢,, » = M,1+ 1,---, M,., are sufficiently 
close together that they may be regarded as equal 
insofar as variation of the elements is concerned. This 
will be expressed 


Vie) Cs a | 
o= Moit+1,--:,M, 


(49) 


Henceforth, the time parameter occurring in the 
argument of an element or element estimate will be 
written 7’, 


It will be supposed that error-free observations are 
given by 


= |e Ete 
Mi1+1,---,M, (50) 


s = 1, 2,.-.-- 


I 


m 


and that the variation of the elements with respect 
to time is given by 


(Ts) = Fi{e(T,), T,, Te). (51) 


Now consider the following stagewise procedure for 
processing the observations (the main reason for 
vhich is its adaptability to cases in which the predic- 
ion functions f“ or F; are pepe oily known): 

The element estimate vector P*(T,) is obtained by 
Minimizing with respect to P* (T,) the quadratic form 


Ni 
= we — gPt(Ty), 4) 
| {PF — gIP*(T,), ol}. (62) 
ee Ts = 1 2. 


b] 


CO Sle RECT eT ar (53) 


Beso IPT ,), Tei T eal (54) 


The quantities P*(T,) will be defined for s > 1 
pelow. To do this, we define sequences of matrices 
B®”, B,°, vy, v4 as follows: 


} Ny 
BY, = De me bP (1), tlbAP*(11), tl (55) 
, B= 
Poa oe ice (56) 
) Ox; b) 
Vo o= (BY, os = 1, 2, (57) 
oe) nie. (fe eee. s| 
k,l=1 (58) 
Gale? Lele Teal 
Gulz,7,, 7) =! f& 7,,7I (59) 
Ox; 
B® =0,- s=1 (60) 
ISG oe ee ee 
M; 
ee = By Dm bs 1P(T.), bb; [P (F2), &) (61) 


b,y= 
Ms—i+1 
Then P*(T,) is obtained for s > 1 by minimizing, 
ith respect to P*(T,), the quadratic form 


fee Set, P.'(7,)} 


4,j=1 


ey (Lj Pi ,)} 


Ms 
+ DS me [FY — g [PT (Ts), tal} 


{F’ — g'|P"(T.), tl} 


Scrutiny of Eqs. (52)—(62) reveals that a stagewise 
kmoothing procedure has been completely defined. 
Now, in order to give the first order error equations 
for the resulting estimates, define the following matrices: 


Ni 
B,i{Ti] = nv bd: [x(T1), bb; {x( 71), tl (63) 


B,y=1 


(where b;’(x, t,) is defined as in (56)); 
¥+17's] =F {By [T.J}~ (64) 
fOrAL SS <P esehan. 


Vil T toy Gulx(T.), Tr, ? T'| 


Te 


cd 
i 


(65) 
-@lj[a(T.), 1s ) TW+nT's) 


(where @,;(z, T; , T’) is defined as in (59)); 
BE ety eer (66) 
and 


Byi{Ts] = BilTs] 


4 Ope e(7), blo/ee(T.), 4). 
TaeOR 
Also define 
eT) = BAT IPS (P.) — 27) 
as (68) 
+ 3 a oFe(T.), bl" — gfe), tb. 
wee 


Then, the first order dependence of P*(T,) — x(T;) 
on the observation errors is: 


PY(T.) — «(Py ={B AP alt (69) 


Also, for T,; < T < Tsi1, the first order dependence 
of P(T) — x(T) is 


PAT) — 2,(T) 
= » @ife(T.), 


Further, suppose that (n,,) represents the inverse 
covariance matrix of {7" — g"[x(t,), tJ}. Then, B,[7,] 
is the inverse covariance matrix of {P;"(7,) — x:(T.)}, 
and B[T] is the inverse covariance matrix of 
{Pil ae De — ae 

In the third section an alternative method was 
described for dealing with time-varying elements. 
That method involved a reduction to the case of con- 
stant elements. 

Choose any time 7, 7. S T° < Tisv, and regard T 
as now fixed. Define a constant element vector x by: 
x = x(T). Let P*(T) represent the estimate of x = 
x(T') obtained by the method described in the third 
section, based on the first 17, observations. 

Then, it can be verified that, to first order, P*(7') = 
P(T), where P(T') is defined by (53). That is, the 
stagewise procedure described in (52)—(62) yields 
an estimate P(7') for x(7') having the same first order 
dependence on F* — g'[x(t,), tJ], » = 1,--:, M., as 
the estimate derived by the method of Section III. 

The stagewise procedure described above actually 
goes through the following steps: P*(T.) represents an 
estimate of x(7’,), based on the first 17, observations; 


T) TYR AG ST) te 
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fort. < f < Diy, PT) is obtained by (53), that 
is, simply by prediction from P*(T,) according to the 
functional relation by which the elements are known 
to vary; P (T's41) represents the estimate Ole haa) 
just before the observations 2”, » = M,+1,:-:,Meu 
are processed. 

The particular method chosen above to define the 
matrices B°—i.e., by (61)—is not the only one 
possible. One could, for example, define matrices 
B® for s > 1 by using (61) with P’(T,) instead of 
P*(T,) in the arguments of b;” and b;’. Then one could 
compute BS”, for s > 1, before computing PAGE: 

This would be convenient for purposes of using 
(20), (23), and (24) for a first order determination of 
P (faa (fone would put p = PT); 
PoP (Ty, Pl, PAL) = (BY yet IPP), 
with p*® defined by the appropriate modification of 
(20) and (23). 

In practical applications, for example, in satellite 
observations where perturbation forces are imper- 
fectly known, the functions ; may not be known ex- 
actly. Furthermore, the inaccuracy in knowledge of 
§{v(T), T, T’} will depend on the time difference 
T’ — T—..e., on how far ahead you are predicting the 
elements. 

Suppose, for example, that the prediction functions 
used in the stagewise procedure are g", 5;, but that 
the correct functions should be h*, 5,*. Let &” still 
represent the observation errors—i.e., &*° = F* — h'— 
and Jet, tor = ¢, , 


é"(T) = LAS en); dis bal, sie a2 
Te id), bal, ta} 
a g{File(T), iN tu, Pia, 
eel) fir bul, ty} 


Suppose that, in the case g* = h', §; = 5,*, the first 
ordermetror im BCT). TT, = Ti <oT's.. could be, ex 
pressed 


(71) 


Ms 
P(T) — a(T) = D) Ter) (72) 
B=1 


ou 
No 
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Then, in the case g" # h", §; # 5.*, one would have 
Ms ; 
PT) —2(T) = » re(T)(e’ + e(T)] (% 
= 


provided 6“(7’) and &* are sufficiently small. 

The errors 6“(7') would, in general, grow as T it 
creases. In such cases, it would be desirable to modif 
the smoothing procedure so as to give more recem 
observations greater weight, and continuously 
diminish the weight given to past observations. There 
are a number of possibilities for accomplishing this 
one way, for instance, would be to define the matrices 
B® in (60)-(62) by 


BP = WO lagrrrs” (008 
where A” and $” are < 1. § 


- 
This would mean that the effective smoothing matrh 


for the observations would not be the matrix (7) 
appearing in (52)—(62); in fact, there would really 
be no fixed smoothing matrix. Also, the first ordei 
error equations (63)—(70) would have to be modifiec 
appropriately. 

There is an analogy between this class of smoothing 
procedures and the process of filtering of signals 
One might regard the smoothing procedures as filters 
with input 2(7') and output P(7T). The case where 
the f" and §; are perfectly known is equivalent to 2 
constant signal input to the filter, in which case the 
filter should have infinite memory to provide maxi. 
mum smoothing of observation errors. If, say, 5; are 
not exactly known, this corresponds to the existence 
of an unpredictable time varying component of the 
filter input; in this case the filter memory must be 
reduced (its ‘bandwidth’ increased); the best ‘time 
constant’ is related to the rate of variation of the un- 
predictable part of x(7’). One could also, in effect 
make the time constant different for the different 
elements. 
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